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Estimating  the  Probability  of  a  Diffusing  Target  Encountering  a 

Stationary  Sensor 

by 

James  N.  Eagle 

Department  of  Operations  Research 

Naval  Postgraduate  School 

Monterey,  CR  93943 

In  this  report  two  expressions  are  given  for  the  probability  of  a 
diffusing  target  avoiding  detection  to  time  t  by  a  sensor  fixed  at  the 
center  of  a  square  region  R.  The  target  is  constrained  to  aluuays  remain 
within  R.  The  first  expression  results  from  an  approximation  to  the 
exact  solution  of  the  diffusion  equation,  and  the  second  from 
experimentation  with  a  Monte  Carlo  simulation  of  the  diffusion  process. 

The  sensor  considered  is  a  definite  range  law  or  "cookie  cutter" 
detector.  For  such  a  sensor,  there  is  a  specified  range,  R,  beyond  which 
detection  is  impossible  and  inside  which  detection  is  certain. 

1.  Background 

This  work  was  begun  with  the  intention  of  developing  for  the  Naval 
war  College,  Newport,  RI  a  simple  expression  for  the  probability  of  a 
diffusing  target  avoiding  detection  by  a  sensor  conducting  a  moving, 
systematic  search  of  the  area  R.   It  was  soon  realized  that  the  special 
case  of  a  stationary  sensor  had  to  be  first  addressed,  and  that  this 
simpler  problem  was  indeed  nontrivial. 

The  importance  of  the  stationary  sensor  problem  goes  beyond  being  a 
limiting  case  of  the  moving  sensor  problem.   It  may  be  used  to  provide 
estimates  of  the  rate  at  which  randomly  moving  targets  will  encounter 


stationary  objects  with  extended  fields  of  influence,  such  as  fixed 
acoustic  sensors,  sonobuoys,  or  possibly  mines. 

This  stationary  sensor  problem  might  appear  deceptively  simple  at 
first  glance,  flfter  all,  it  could  be  argued,  this  problem  is  equivalent  to 
that  of  a  searcher  with  a  cookie  cutter  sensor  of  range  R  conducting  a 
random  search  for  a  stationary  target,  find  Koopman[1946]  and  [1980] 
argued  that  the  probability  of  nondetection  to  time  t  is  exp(-2Rvt/fl), 
where  v  is  the  speed  of  the  random  search.  The  problem  is,  of  course, 
what  to  use  for  v  when  the  searcher's  path  is  a  diffusion.  One  of  the 
results  of  this  work  is  an  expression  for  the  equivalent  speed  of  such  a 
"diffusion  search"  of  a  square  area. 

The  initial,  experimental  results  for  the  problem  here  addressed  were 
obtained  by  Sislioglu[1981].  Sislioglu  conducted  Monte  Carlo  simulations 
with  different  target  diffusion  constants  D,  d\red  sizes  R,  and  sensor 
detection  ranges  R.  He  observed  that  when  the  initial  distribution  of  the 
target  was  uniform  over  R,  the  probability  of  nondetection  to  time  t, 
PND(t),  was  given  approximately  by 

(1  -  TTR2/R)  exp(-247  RDt/R15).  (1) 

In  this  report,  some  analytical  support  for  Sislioglu's  results  is  offered. 
Rlso  a  slight  modification  of  (1}  is  suggested  which  agrees  more  closely 
with  theory  and  expermental  results  when  the  area  of  the  detection 
disk  approaches  the  ar^d  of  the  search  region  R.  , 

2.  The  Diffusion  Equation 

The  probability  density  p(t)  of  a  particle  undergoing  diffusion  in  any 
coordinate  system  must  satisfy  the  diffusion  equation 

(D/2)  V2p  =  3p/9t,  (2) 


where  D  is  the  diffusion  constant,  and  V2  is  the  Laplacian  operator.   In 
Cartesian  and  polar  coordinates,  respectively,  (2)  becomes 

ID/2)  02p/3x2  +  32p/3y2)  =  3p/at,  and  (3) 

(d/2)  o2p/ar2  +  (i/r2)02p/ae2)  +  d/r)Op/ar})  =  ap/at.  H) 

To  find  a  unique  function  p(x,y,t)  or  p(r,8,t)  satisfying  these  partial 
differential  equations,  spatial  and  temporal  boundary  conditions  must  be 
specified.  Defining  R  to  be  a  square  with  sides  of  length  L  in  the  first 
quadrant,  the  boundary  conditions  in  Cartesian  coordinates  d^e 

9p/ax|x=o  =    9p/9*|x=L  =  o,  (5) 

3p/ay|y=0  -    9P/9y|y=L  -  o,  (6) 

p(x,y,t)  =  0,  ((x  -  L/2)2+(y  -  L/2)2  )1/2  <  R,  and  (?) 

p(x,y,0)  =  1/R,  (lx  -  L/2)2+(y  -  L/2)2  I1'2  >  R.  (8) 

Equations  (5)  and  (6)  ensure  that  none  of  the  target's  probability  mass 
"escapes"  from  R.  That  is,  the  boundaries  of  R  are  reflecting.  Equation 
(7)  requires  that  p(x,y,t)  be  0  on  the  detection  disk.  Rnd  (8)  ensures  that 
the  initial  distribution  of  the  target  over  the  search  region  is  uniform  and 
integrates  to  (R-ttR2)/R  (the  probability  that  the  target  is  not  detected 
at  time  0). 

For  any  particular  instance  of  the  problem,  finding  a  p(x,y,t) 
satisfying  (3)  and  the  boundary  conditions  is  not  difficult  using  numerical 
methods.    Such  procedures  are  routinely  used  in  heat  transfer  problems 
to  solve  the  diffusion  equation  (called  the  conduction  equation  or  the 
Fourrier  equation  in  the  physics  and  engineering  literature)  to  determine 
the  temperature  distribution  across  imperfectly  conducting  solids.   In 
fatt,  Pitts  and  Sissom(1977]  give  the  example  of  a  heated  pipe  in  a 
square  block  of  insulating  material  as  one  where  the  isotherms  can  be 
accurately  estimated  by  hand  plotting. 


RLthough  the  problem  is  not  hard  to  solve  numerically,  the  square 
boundary  of  R  combined  with  a  circular  sensor  tend  to  make  the 
analytical  solution  difficult  to  obtain,  find  without  an  analytical 
solution,  it  may  be  impossible  to  establish  Sislioglu's  observations  in 
general.    Making  a  change  to  the  geometry,  however,  allows  an 
analytical  solution.  Specifically,  if  the  search  region  R  is  assumed  to  be 
a  disk  of  area  R  instead  of  a  square,  then  an  exact  solution  in  polar 
coordinates  is  possible. 

It  is  noted  that  the  ray  solution  method  described  by  Mangel  [1981  ] 
could  presumably  be  used  to  solve  the  diffusion  equation,  at  least 
approximately,  for  a  square  search  area.  Such  a  solution  was  not 
attempted  since  an  exact  solution  was  available  for  the  circular  search 
area  case. 

3.  The  Solution 

The  disk-within-a-disk  geometry  has  a  radial  symmetry,  thus  reducing 
the  problem  to  one  dimension,  r.  The  new  problem  is  to  find  a  function 
p(r,t]  satisfying 

(D/2)  (d2p/9r2  +  (1/r)Op/3r))  =  9p/9t, 

subject  to 

ap/ar|r=Rfl  =  o, 

p(r,t)  =  0,  r«R, 

p(r,0)  =  1/fl,  r>R,  " 
where  Rfl  =  (R/ir)1/2is  the  radius  of  the  transformed  (i.e.,  circular)  search 
area  R.  This  problem  has  been  solved  in  the  physics  literature. 


Muskat[1937]  and  Muskat[1934]  (a  paper  investigating  the  production 
rate  of  oil  wells)  give  the  solution  as 

00 

p(r,t)    =  -lff/fl)  £  kn  U(ornr)  exp(-D  <*r2 1/2)  (9) 

n=1 

where 

kn  =  J0(ornR)  J^nPn)  /  U02l<*nR)  -  Ji2(ornRfi)),  (10) 

U(anr)  =  Y^nRfl)  J0(arrir)  -  J1(»r,Rfi)  Yrfar„r),  (11) 

and  an  is  the  n*  positive  root  of  U(a-R).  (That  is,  the  nth  smallest 
positive  value  of  a  satisfying  U(orR)  =  0.)  RLso,  J0,  J*  Y0,  and  Yj  &re 
respectively  Bessel  functions  of  the  first  kind  order  0,  first  kind  order  1, 
second  kind  order  0,  and  second  kind  order  1. 

The  solution  (9H11)  does  not  appear  particularly  easy  to  evaluate  or 
interpret,  being  an  infinite  sum  of  Bessel  functions.  However  for  large  t 
the  solution  simplifies  considerably.  Rs  t  becomes  large,  (9)  becomes 
l-irk,  /R)  UU*,r]  exp(-D  c*,2t/2],  (12) 

iut\ere  or1  is  the  smallest  positive  value  of  &  such  that  UlctR)  =  0.  The 
other  exponential  terms  are  ignored  since  they  involve  larger  roots  of 
UlceR].  This  means  that  as  t  becomes  large,  the  decrease  in  p(r,t)  for 
constant  r  becomes  exponential  at  a  rate  of  D<*  2/2. 

When  p(r,t)  is  specified,  PND(t)  is  then  given  by 

I    I  p(r,t)rdrd6.  (13) 

0     R 


So  for  Large  t,  PND(t)  becomes 

-  t2tf2k1/R)  {J  U{<*,r)  r  dr}  expl-D  c^t/2)  (14) 

R 

=      -  (2lT%<|  R/CRori  DiJ^tt!  Rfll  Yjt*,  R)  -  Y^RJ  J^  R)}  exp(-D  <*,2t/2),      (15) 
K    exp  (-D  *,2t/2), 
Luhere  K  is  a  constant  tuhich  depends  on  the  problem  geometry. 
Evaluation  of  the  integral  in  (14)  is  straightforward  given  the  change  of 
variable  u= o>^1  r  and  the  identities 

|x  J0(x)  dx  =  x  Jjlx)  and     Jx  YCl(x)  dx  =  x  Y^x). 
Thus  for  large  t,  PND(t)  decreases  exponentially  at  the  same  rate  as 
p(r,t).  So  Sislioglu's  observation  of  exponentially  decreasing  PND(t)  (or, 
equivalently,  a  constant  detection  rate)  appears  asymptotically  correct 
for  large  enough  t. 

Using  Muskat's  dimensionless  notation,  we  can  simplify  the  solution 
somewhat  by  defining 

x,=  <*tR  and  p  =  Rft/R. 
Then  for  Large  t,  PND(t)  can  be  written  as 

KexpKix^t/^R2)),  (16) 

where  x,  is  the  first  positive  x  satisfying 

Y^xpJJofxJ-J^xpjYofx)  =  0;  (1?) 

K  is  given  by 

-(271k,  /(xlP2))  {J^p)  YJlxJ  -  Yrfx,p)  JM};  (18) 

and  k,  is 

JrfxJ  Jjx,  p)  /  (J02(x,)  -  jflx,  p)J.  (19) 


In  the  next  section,  experimentally  determined  PNDl't)  from  a  Monte 
Carlo  simulation  of  the  diffusion  process  will  be  compared  with  (16H19), 
Sislioglu's  approximation  (1),  and  a  slight  modification  of  (1). 


4.  Simulation  Results  and  Conclusions 

Figure  1.  is  a  plot  of  PND(t)  determined  by  Monte  Carlo  simulation  of 
the  diffusion  process  for  a  square  area  R  of  size  10,000  square  nautical 
miles  (nm2),  a  diffusion  constant  D  of  100  nm2/hour,  and  detection  radii  R 
varying  from  28.21  nm  (p  =  2)  to  3.76  nm  (p=15).  The  time  increment  for 
these  simulations  was  1  minute,  which  resulted  in  the  x  and  y 
displacements  of  the  target  in  each  increment  being  selected  from 
independent  normal  distributions  of  mean  0  and  variance  100/60. 
Figure  1.  illustrates  the  degree  to  which  the  decrease  in  PND(t)  is 
exponential.  Plotted  on  a  log  scale,  PND(t)  appears  very  nearly  linear. 
For  small  t,  however,  the  decrease  in  PND(t)  is  faster  than  exponential. 
This  is  seen  more  clearly  in  Figure  2,  which  is  an  enlargement  of  the 
upper  left-hand  corner  of  Figure  1. 

Sislioglu  reported  that  PND(t)  can  be  approximated  by 

(1  -  TTR2/R)  exp(-  24.7  RDt/fl  15L  (20) 

The  simulation  results  reported  here  indicate  that  a  somewhat  better 
approximation  when  p  approaches  1  is . 

(1  -  7TR2/R)  exp(-  24.7  R0t/(R  -  TTR2)15).  (21) 

That  is,  R  in  (20)  is  replaced  with  R  -  7TR2    Figure  3.  compares  PNO(t) 
calculated  by  simulation  with  the  estimates  given  by  (21),  (20),  and  the 
asymptotic  estimate  (16H19).  The  simulation  data  in  Figure  3.  are  the 
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same  as  in  Figures  1.  and  2.  except  that  only  results  for  p  of  2,  3,  and  6 
are  shown.   In  Figure  3.  the  simulated  PND(t)  is  shown  as  a  solid  line, 
while  the  estimated  PND(t)  is  shown  as  dashed. 

Figure  3.(0  indicates  that,  especially  for  small  p,  PfND(t)  determined 
by  simulation  does  not  decrease  as  rapidly  as  predicted  by  (16H19).  The 
explanation  is  believed  to  be  that  for  p  sufficiently  close  to  1,  a  circular 
search  region  does  not  provide  a  reasonable  approximation  of  a  square 
region  of  the  same  area.  To  test  this  explanation,  the  three  simulations 
plotted  in  Figure  3.  were  repeated  with  circular  search  areas.  The 
results,  shown  in  Figure  3.(D).,  indicate  a  closer  agreement  between  the 
theoretical  and  observed  data.  But  still  the  fit  is  not  exact.  This  is 
somewhat  disturbing,  but  might  be  explained  by  the  discrete  manner  in 
which  the  diffusion  path  is  simulated.  The  simulated  diffusion  path  is 
approximated  by  a  series  of  points.  Detection  must  occur  exactly  at  the 
points  and  can  not  occur  between  them.  The  distance  between  adjacent 
points  is  the  random  variable  (AX2+AY2)1/2,  where  AX  and  AY  are  the 
independent,  normally  distributed  x  and  y  displacements.   It  is  possible 
for  the  simulated  path  to  jump  across  the  edge  of  the  detection  disk 
without  achieving  detection,  even  though  the  line  segment  connecting 
the  points  lies  partly  on  the  disk.  This  will  tend  to  reduce  the  simulated 
detection  rate  below  that  of  the  diffusion  being  approximated. 

Figure  4.  shows  a  plot  of  Xjvs.  p  for  values  of  p  from  2  to  15.  By 
using  (16),  these  values  of  y^  determine  the  theoretical  asymptotic 
detection  rate  as 

D  x12/(2R2).  (22) 

Table  1.  lists  the  asympolic  detection  rates  determined  by  (22),  (21),  (20), 


or 


-   ao    C£ 


O 

X 

q: 


and  an  overall  (i.e.,  from  time  0]  best  fit  rate  calculated  by  least- 
squares  fitting  of  the  simulation  data. 


xjp)  Dxt2/(2R2]     24.7  RD/(FMTR2}1-5     24.7  RD/fl15      Simulation 

(Eq.  22]  (Eq.  21)  (Eq.  20)  Best  Fit 


2 

1.361 

.116 

.107 

.0697 

.102 

2.5 

.866 

.0736 

.0724 

.0557 

.0665 

3 

.626 

.0554 

.0554 

.0465 

.0507 

4 

.389 

.0381 

.0384 

.0384 

.0388 

6 

.218 

.0269 

.0242 

.0232 

.0250 

8 

.147 

.0217 

.0178 

.0174 

.0202 

10 

.111 

.0194 

.0141 

.0139 

.01704 

15 

.0661 

.0154 

.00935 

.00929 

.0137 

Table  1.  Detection  Rates  for  Various  p. 

The  simulation  data  suggest  the  following  two  conclusions: 

a.  Equations  (16H19)  give  a  reasonable  estimate  PND(t),  but  the  fit 
deteriorates  as  p  decreases  to  1. 

b.  For  small  p,  (21)  gives  a  better  fit  to  the  data  than  does  (20). 
For  large  p,  both  (20)  and  (21)  underestimate  the  observed  asymptotic 
detection  rate. 

LUe  conclude  with  a  few  comments  on  random  search.  Rs  mentioned 
earlier,  Koopman's  random  search  model  predicts  a  detection  rate  of 
2Rv/R  for  searcher  with  speed  v  and  detection  range  R  conducting  a 
random  search  of  an  area  of  size  R.   It  seems  reasonable  that  a 
diffusion  path  should  be  "random"  enough  for  this  model  to  be 
appropriate.   In  fact,  we  have  seen  that  the  detection  rate,  while  not 


constant,  approaches  a  constant  value  for  Large  t.  Setting  2Rv/R  equal 
to  the  exponential  term  in  (21)  and  solving  for  speed  v  gives 

12.35  RD/IR  -  ffR2)1-5 
as  an  approximate  equivalent  speed  for  a  searcher  conducting  a 
"diffusion  search"  of  a  square  area. 
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